5 Hazard Flags and BER

1 Load Data

load(file='./source/in_vitro/tier1_NAMs_apcra_pro_02262024.RData')
load(file='./source/chem/apcra_chem_ad.RData')
asnm.tier1.all.invivo <- read.xlsx('./output/asnm_tier1_all_invivo.xlsx', colNames = TRUE, sheet=1) %>% as.data.table()
seem3 <- fread('./source/exposure/SupTable-all.chem.preds-2018-11-28.txt')

# I tried using this RData from Git but it did not contain any of our DTXSIDs:
#load('./source/exposure/new.chem.preds-2019-12-13.RData')
#seem3 <- chem.preds

colnames(mega.mc5) # this is going to have the hazard information we need

##  [1] "spid"                "chid"                "casn"               
##  [4] "chnm"                "dsstox_substance_id" "code"               
##  [7] "aeid"                "aenm"                "m5id"               
## [10] "m4id"                "bmad"                "resp_max"           
## [13] "resp_min"            "max_mean"            "max_mean_conc"      
## [16] "max_med"             "max_med_conc"        "logc_max"           
## [19] "logc_min"            "cnst"                "hill"               
## [22] "hcov"                "gnls"                "gcov"               
## [25] "cnst_er"             "cnst_aic"            "cnst_rmse"          
## [28] "cnst_prob"           "hill_tp"             "hill_tp_sd"         
## [31] "hill_ga"             "hill_ga_sd"          "hill_gw"            
## [34] "hill_gw_sd"          "hill_er"             "hill_er_sd"         
## [37] "hill_aic"            "hill_rmse"           "hill_prob"          
## [40] "gnls_tp"             "gnls_tp_sd"          "gnls_ga"            
## [43] "gnls_ga_sd"          "gnls_gw"             "gnls_gw_sd"         
## [46] "gnls_la"             "gnls_la_sd"          "gnls_lw"            
## [49] "gnls_lw_sd"          "gnls_er"             "gnls_er_sd"         
## [52] "gnls_aic"            "gnls_rmse"           "gnls_prob"          
## [55] "nconc"               "npts"                "nrep"               
## [58] "nmed_gtbl"           "hitc"                "modl"               
## [61] "fitc"                "coff"                "actp"               
## [64] "modl_er"             "modl_tp"             "modl_ga"            
## [67] "modl_gw"             "modl_la"             "modl_lw"            
## [70] "modl_rmse"           "modl_prob"           "modl_acc"           
## [73] "modl_acb"            "modl_ac10"           "resp_unit"          
## [76] "conc_unit"           "mc6_flags"           "flag.length"        
## [79] "use.me"              "ac50_uM"             "acc_uM"             
## [82] "asnm"

2 Calculating BERs

2.1 Examine and add SEEM3

• Here are the 8 APCRA substances for which there are no SEEM3 data reported.
• They appear generally to be salts, with a few that are not.

seem.apcra <- seem3[dsstox_substance_id %in% apcra.list[,DTXSID]] # only 193...

apcra.list[!DTXSID %in% seem.apcra$dsstox_substance_id]

• Fifteen chemicals had missing SEEM3 U95 values, see below.

seem.apcra[is.na(seem3.u95) & !is.na(seem3)]

• For these, the median SEEM3 estimate was substituted.
• log10-SEEM3 values were also calculated for log10 BER computation.

seem.apcra[is.na(seem3.u95)& !is.na(seem3), seem3.u95 := seem3]
seem.apcra[,seem3.log10 := log10(seem3)]
seem.apcra[,seem3.u95.log10 := log10(seem3.u95)]

• Adding in SEEM3 to the sheet and calculating BERs.

asnm.tier1.all.invivo$seem3.u95.log10 <- seem.apcra$seem3.u95.log10[match(asnm.tier1.all.invivo$dsstox_substance_id,
                                                seem.apcra$dsstox_substance_id)]

asnm.tier1.all.invivo$seem3.log10 <- seem.apcra$seem3.log10[match(asnm.tier1.all.invivo$dsstox_substance_id,
                                                seem.apcra$dsstox_substance_id)]

asnm.tier1.all.invivo[,ber.targeted := median(aed50.atg,
                                    aed50.bsk,
                                    aed50.ccte,
                                    aed50.nvs, na.rm=TRUE) - seem3.u95.log10]


asnm.tier1.all.invivo[,ber.httr := median(aed50.httr.heparg,
                                       aed50.httr.u2os,
                                       aed50.httr.mcf7, na.rm=TRUE) - seem3.u95.log10]


asnm.tier1.all.invivo[,ber.htpp := ifelse(!is.na(aed50.htpp.u2os), aed50.htpp.u2os - seem3.u95.log10, NA)]


asnm.tier1.all.invivo[,ber.astar := median(aed50.astar.beas2b,
                                    aed50.astar.hek293,
                                    aed50.astar.hepg2,
                                    na.rm=TRUE) - seem3.u95.log10]

asnm.tier1.all.invivo[,ber.med.aed50 := med.aed50 - seem3.u95.log10]

review <- 
asnm.tier1.all.invivo[ber.med.aed50 < 4]

• In general, BER was NA when SEEM3 was missing, except for one silane with no median AED50.

asnm.tier1.all.invivo[is.na(ber.med.aed50), c('chnm','med.aed50','seem3.u95.log10')] # these are missing exposure values and/or bioactivity values

• 56 APCRA pro only chemicals with BER <4
• 40 APCRA pro only chemicals with BER <3

asnm.tier1.all.invivo <- asnm.tier1.all.invivo[order(-ber.med.aed50)]
asnm.tier1.all.invivo[ber.med.aed50 < 4 & apcra.pro.only==1, 
                      c('chnm','apcra.pro.only','med.aed50','seem3.u95.log10','ber.med.aed50')] #56 rows

• As a group the silanes generally have caution on their AQC.
• Presence of silanes in the low BER group is therefore maybe associated with more uncertainty.

asnm.tier1.all.invivo[grep('silane',chnm),c('chnm','dsstox_substance_id','T0','T4','Call','aqc_iv_pass','aqc_indicator')]

2.2 Prepare BER plotting

asnm.tier1.all.plot <- melt.data.table(asnm.tier1.all.invivo,
                                id.vars = c('dsstox_substance_id',
                                            'chnm',
                                            'CASRN',
                                            'apcra.pro.only',
                                            "ber.targeted", 
                                            "ber.httr", 
                                            "ber.htpp",
                                            "ber.astar",
                                            "ber.med.aed50",
                                            "aqc_indicator"),
                                measure.vars = c("log.repdose.any.5th",
                                                 "aed50.atg", 
                                                 "aed50.bsk", 
                                                 "aed50.ccte", 
                                                 "aed50.nvs", 
                                                 "aed50.stm", 
                                                 "aed50.httr.mcf7", 
                                                 "aed50.httr.u2os", 
                                                 "aed50.httr.heparg", 
                                                 "aed50.htpp.u2os", 
                                                 "aed50.astar.beas2b", 
                                                 "aed50.astar.hepg2", 
                                                 "aed50.astar.hek293", 
                                                 "med.aed50", 
                                                 "min.aed50", 
                                                 "seem3.u95.log10", 
                                                 "seem3.log10"),
                                variable.name = 'types',
                                value.name = 'values')

asnm.tier1.all.plot <- asnm.tier1.all.plot[order(-ber.med.aed50)]

ber.all <- ggplot(data=asnm.tier1.all.plot[!is.na(ber.med.aed50)], 
            aes(x=reorder(factor(chnm),-ber.med.aed50), 
                y=values,
                group = factor(types),
                shape = factor(types), 
                color = factor(types)))+
  geom_point(size=3)+
  scale_colour_viridis(discrete=TRUE,
                       name = 'Values',
                      breaks=c("log.repdose.any.5th", "aed50.atg", 
"aed50.bsk", "aed50.ccte", "aed50.nvs", "aed50.stm", "aed50.httr.mcf7", 
"aed50.httr.u2os", "aed50.httr.heparg", "aed50.htpp.u2os", "aed50.astar.beas2b", 
"aed50.astar.hepg2", "aed50.astar.hek293", "seem3.u95.log10", "seem3.log10"),
                      labels=c('5th%-ile POD All',
                               'ATG AED50',
                               'BSK AED50',
                               'MEA AED50',
                               'NVS AED50',
                               'STM AED50',
                               'HTTr MCF7 AED50',
                               'HTTr U2OS AED50',
                               'HTTr HepaRG AED50',
                               'HTPP U2OS AED50',
                               'ASTAR BEAS2B AED50',
                               'ASTAR HepG2 AED50',
                               'ASTAR HEK293 AED50',
                               'SEEM3 U95',
                               'SEEM MED'
                            )) +
  scale_shape_manual(name = 'Values',
                     breaks=c("log.repdose.any.5th", "aed50.atg", 
"aed50.bsk", "aed50.ccte", "aed50.nvs", "aed50.stm", "aed50.httr.mcf7", 
"aed50.httr.u2os", "aed50.httr.heparg", "aed50.htpp.u2os", "aed50.astar.beas2b", 
"aed50.astar.hepg2", "aed50.astar.hek293", "seem3.u95.log10", "seem3.log10"),
                     values=c(15,16,17,18,8,9,15,16,17,18,8,9,15,16,17,18,8,9),
                     labels=c('5th%-ile POD All',
                               'ATG AED50',
                               'BSK AED50',
                               'MEA AED50',
                               'NVS AED50',
                              'STM AED50',
                               'HTTr MCF7 AED50',
                               'HTTr U2OS AED50',
                               'HTTr HepaRG AED50',
                               'HTPP U2OS AED50',
                               'ASTAR BEAS2B AED50',
                               'ASTAR HepG2 AED50',
                               'ASTAR HEK293 AED50',
                               'SEEM3 U95',
                               'SEEM MED'))+
  xlab('')+
  #ylab('log10-mg/kg/day')+
  theme_bw() +
  theme(
    #axis.text.x = element_text(angle=90, vjust=0.5),
        axis.title.y = element_text(size=12, face='bold'),
        axis.title.x = element_blank(),
        axis.text.y = element_blank(),
        legend.position = 'bottom')+
  scale_y_continuous(breaks=seq(-13,5,2))+
  guides(colour=guide_legend(nrow=5))+
      annotate("rect", 
           ymin = -9, 
           ymax = 7, 
           xmin = 134, 
           xmax = 190, 
           fill = "red", 
           alpha = 0.1, 
           size = 1, 
           color = "red")+
    annotate("text",
           x = 188,
           y=-11,
           label="A",
           size=10,
           hjust=0,
           vjust=1)+
  
  coord_flip()
#+
#  theme(axis.text.y = element_text(colour=retro))


ber.all

• Reduce complexity

asnm.tier1.b <- asnm.tier1.all.plot[ber.med.aed50 <4 & aqc_indicator==1 & apcra.pro.only==1]
#unique(asnm.tier1.b$dsstox_substance_id) #13 substances
#retro <- ifelse(asnm.tier1.a$apcra.ret==1, "red","black")
#retro.face <- ifelse(asnm.tier1.a$apcra.ret==1, "italic","plain")

ber.simple <- ggplot(data=asnm.tier1.b[ types %in% c("log.repdose.any.5th", 
                                                                                              "med.aed50",
                                                                                              "seem3.log10",
                                                                                              "seem3.u95.log10")], 
            aes(x=reorder(factor(chnm),-ber.med.aed50), 
                y=values,
                group = factor(types),
                shape = factor(types), 
                color = factor(types)))+
  geom_point(size=3)+
  scale_colour_manual(values=c('#481567FF',
                               '#2D708EFF',
                               '#3CBB75FF',
                               '#95D840FF'),
                       name = 'Values',
                      breaks=c("log.repdose.any.5th", 
                               "med.aed50",
                               "seem3.u95.log10", 
                               "seem3.log10"),
                      labels=c('5th%-ile POD All',
                               "Med AED50",
                               "SEEM3 U95",
                               "SEEM3"
                            )) +
  scale_shape_manual(name = 'Values',
                     breaks=c("log.repdose.any.pod", 
                               "med.aed50",
                               "seem3.u95.log10", 
                               "seem3.log10"),
                      labels=c('5th%-ile POD All',
                               "Med AED50",
                               "SEEM3 U95",
                               "SEEM3"),
                     values=c(15,16,17,18))+
  xlab('')+
  #ylab('log10-mg/kg/day')+
  theme_bw() +
  theme(
    #axis.text.x = element_text(angle=90, vjust=0.5),
        axis.text = element_text(size=10),
        axis.title.y=element_text(size=14,face='bold'),
        axis.title.x = element_blank(),
        legend.position = 'bottom')+
  scale_y_continuous(breaks=seq(-9,5,2))+
  guides(colour=guide_legend(nrow=2))+
  annotate("text",
           x = 43,
           y=-8,
           label="C",
           size=10,
           hjust=0,
           vjust=1)+
  coord_flip()
  #theme(axis.text.y = element_text(colour=retro, face=retro.face))

  
ber.simple

bplusc <- plot_grid(ber.complex, ber.simple, rel_heights = c(1,1), rel_widths = c(1,1), nrow=2)

plot_grid(ber.all, bplusc, rel_widths = c(1,1), rel_heights = c(1,1), ncol=2)

layout <- matrix(c(1,1,2,2,2,
                   1,1,2,2,2,
                   1,1,3,3,3,
                   1,1,3,3,3), nrow=4, byrow=TRUE)

# multiplot was obtained from
# http://www.cookbook-r.com/Graphs/Multiple_graphs_on_one_page_%28ggplot2%29/
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  require(grid)
  
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  
  numPlots = length(plots)
  
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  
  if (numPlots==1) {
    print(plots[[1]])
    
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

file.dir <- paste("./output", sep="")
file.name <- paste("/Figure_BER_", Sys.Date(), ".png", sep="")
file.path <- paste(file.dir, file.name, sep="")
dir.create(path=file.dir, showWarnings = FALSE, recursive = TRUE)
png(file.path, width = 10000, height = 9000, res=600)
multiplot(ber.all, ber.complex, ber.simple, layout = layout)
dev.off()

## png 
##   2

Table1redo <- asnm.tier1.all.invivo[ber.med.aed50< 4 & apcra.pro.only==1 & is.na(log.repdose.any.25th) & aqc_indicator==1, c('dsstox_substance_id','preferred_name','ber.med.aed50', 'log.repdose.any.25th','med.aed50', 'seem3.u95.log10', 'seem3.log10')]
Table1redo <- Table1redo %>% mutate_if(is.numeric, ~round(., 2))
Table1redo

Table2redo <- asnm.tier1.all.invivo[ber.med.aed50 > 3 & med.aed50 > 2 & pod.ratio > -0.5 & aqc_indicator==1, c('dsstox_substance_id','preferred_name','log.repdose.any.5th', 'med.aed50', 'ber.med.aed50', 'pod.ratio', 'aqc_indicator')]

Table2redo <- Table2redo %>% mutate_if(is.numeric, ~round(., 2))
Table2redo

redo_tables <- list('Table1_HighPriorityChems' = as.data.frame(Table1redo),
                    'Table2_LowPriorityChems' = as.data.frame(Table2redo))


write.xlsx(redo_tables, './output/Draft_Table1_Table2_APCRA_Pro_27May2024.xlsx')

2.3 Compare HepaRG POD to other POD

heparg.comp <- asnm.tier1.all.invivo[,c('dsstox_substance_id',
                                        'aed50.httr.heparg',
                                        'aed50.httr.u2os',
                                        'aed50.httr.mcf7',
                                        'med.aed50',
                                        'min.aed50',
                                        'p5.toxval.numeric',
                                        'p25.toxval.numeric'
                                        )]
heparg.comp[, httr.min.aed50.nonheparg := min(aed50.httr.heparg,aed50.httr.mcf7, na.rm = TRUE), by=c('dsstox_substance_id')]
heparg.comp[, httr.heparg.diff := ifelse(!is.na(httr.min.aed50.nonheparg), httr.min.aed50.nonheparg - aed50.httr.heparg, NA), by=c('dsstox_substance_id')]

range(heparg.comp$httr.heparg.diff, na.rm=TRUE) #-3.2 to 0

## [1] -3.199864  0.000000

hist(heparg.comp$httr.heparg.diff, na.rm=TRUE)

heparg.comp.long <- melt.data.table(heparg.comp,
                                    id.vars = c('dsstox_substance_id',
                                                'aed50.httr.heparg'),
                                    measure.vars = c('httr.min.aed50.nonheparg',
                                                     'med.aed50',
                                                     'p5.toxval.numeric',
                                                     'p25.toxval.numeric'),
                                    variable.name = c('POD'),
                                    value.name = c('value'))

heparg.httr <- ggplot(data=heparg.comp.long, 
                 aes(x=aed50.httr.heparg, 
                     y=value))+
  geom_point(aes(x=aed50.httr.heparg, 
                 y=value), alpha=0.5, size=2)+
  coord_cartesian(xlim=c(-4,6), ylim=c(-4,6))+
  theme_bw()+
  geom_abline(color='black')+
  geom_abline(intercept = 0.5, slope = 1, color="dark gray", 
                 linetype="dashed", size=0.75)+
  geom_abline(intercept = -0.5, slope = 1, color="dark gray", 
                 linetype="dashed", size=0.75)+
  xlab('HepaRG HTTr AED50, log10-mg/kg/day')+
  ylab('POD value, log10-mg/kg/day')+
   theme(axis.text = element_text(size=12),
        axis.title = element_text(size=14))+
  facet_wrap(~POD, scales='fixed')+
  theme(strip.background = element_rect(fill='white',size=1),
        strip.text = element_text(size=10, color='black',face='bold'))

heparg.httr

file.dir <- paste("./output", sep="")
file.name <- paste("/Figure_HepaRG_", Sys.Date(), ".png", sep="")
file.path <- paste(file.dir, file.name, sep="")
dir.create(path=file.dir, showWarnings = FALSE, recursive = TRUE)
png(file.path, width = 4000, height = 4000, res=600)
heparg.httr
dev.off()

2.4 SEEM Credible Interval Size vs BER

ber.v.seem <- asnm.tier1.all.invivo[,c('dsstox_substance_id','ber.med.aed50', 'seem3.u95.log10','seem3.log10', 'med.aed50')]
ber.v.seem[ ,cred.int.size := seem3.u95.log10 - seem3.log10]


ber_seem <- ggplot(data=ber.v.seem, 
                 aes(x=cred.int.size, 
                     y=ber.med.aed50))+
  geom_point(aes(x=cred.int.size, 
                     y=ber.med.aed50), alpha=0.5, size=2)+
  #coord_cartesian(xlim=c(-4,6), ylim=c(-4,6))+
  theme_bw()+
  geom_smooth(method='lm',se=FALSE,color='blue',formula=y~x)+
  stat_poly_eq(formula=y ~x,
               eq.with.lhs = "y~`=`~",
               aes(label=paste0("atop(", ..eq.label.., ",", ..rr.label.., ")")),
               label.x.npc = 0.8,
               label.y.npc = 0.8,
               parse=TRUE)+
  xlab('SEEM3 U95 - SEEM3 Median, log10-mg/kg/day')+
  ylab('BER')+
   theme(axis.text = element_text(size=12),
        axis.title = element_text(size=14))
ber_seem

file.dir <- paste("./output", sep="")
file.name <- paste("/Figure_BER_v_SEEMCredInt_", Sys.Date(), ".png", sep="")
file.path <- paste(file.dir, file.name, sep="")
dir.create(path=file.dir, showWarnings = FALSE, recursive = TRUE)
png(file.path, width = 4000, height = 4000, res=600)
ber_seem
dev.off()

## png 
##   2

3 Developmental and reproductive toxicity (DART)

Dev tox flag will be added to the ER/AR flags as a set for DART prediction flag. The dev tox flag consists of Stemina devTox quick-Predict data and an in silico flag from TEST.

3.1 Reference chemicals for developmental toxicity

Using Zurlinden et al 2019 for reference chemicals. Can see that the devtox TEST prediction is not necessarily accurate.

devref <- fread('./source/in_vitro/devtox_reference_zurlinden_2019.csv')
devref <- devref[-26,]
devref[1:17, ref:= 'TP']
devref[18:26, ref:= 'FN']
devref[27:42, ref:= 'TN']
devref[,c('DTXSID','PREFERRED_NAME','ref','DEVTOX_TEST_PRED')]

What is the overlap with APCRA 201 substances? Seemingly only one true positive…and not one captured by Stemina (PTU).

devref[DTXSID %in% apcra.list$DTXSID]

Load mc5 for some positive reference chemicals for DART overall (dev, ER, AR)

dart.ref <- c(
'DTXSID5022308', #genistein
'DTXSID3020465', #diethylstilbestrol
'DTXSID7020182', #bisphenol a
'DTXSID7032004', #flutamide
'DTXSID4022361', #vinclozolin
'DTXSID7021239', #retinoic acid
'DTXSID9022524', #thalidomide
'DTXSID2020634', #5-fluorouracil
'DTXSID6025438', #hydroxyurea
'DTXSID1020194'  #boric acid
)

Here we pull the DART reference chemical mc5 for Stemina from invitrodb v3.5.

dart.mc5 <- tcplPrepOtpt(tcplLoadData(lvl=5, type='mc',fld='aeid', val=c(1691,1858)))
dart.mc5 <- dart.mc5[dsstox_substance_id %in% dart.ref]

3.2 Add dev datasets

• 61 APCRA substances were screened in multi-conc in Stemina (STM).
• An additional 5 substances appeared positive in the single-conc screen but were not re-screened; as such a positive is inferred at the concentration screened in single conc (100 uM).

stm <- mega.mc5[asnm=='STM']
#length(unique(stm$dsstox_substance_id)) #66
stm.mc.dtxsids <- unique(stm$dsstox_substance_id) # 6 of these were added as positives from sc that never advanced to mc (from previous processing)

#stm[hitc==1 & is.na(modl_tp)]
#added.dtxsid <- stm[is.na(modl_prob), dsstox_substance_id]
#stm[dsstox_substance_id %in% added.dtxsid, modl_acc := 2]
#stm[dsstox_substance_id %in% added.dtxsid, acc_uM := 100]

• Reshape the multi-conc data.

# cast data to combine them.

mc5.stm.wide <- dcast.data.table(stm[aeid %in% c(1691,1858)], dsstox_substance_id + chnm ~ aenm,
                                  value.var = 'modl_acc',
                                  fun.aggregate= min)

mc5.dartref.wide <- dcast.data.table(dart.mc5[aeid %in% c(1691,1858)], dsstox_substance_id + chnm ~ aenm,
                                  value.var = 'modl_acc',
                                  fun.aggregate= min)

mc5.stm.wide <- rbind(mc5.stm.wide,
                      mc5.dartref.wide)

• Negatives need to be inferred for the single conc screened substances that were negative.
• 133 chemicals with single conc screening data from the APCRA list.

add.wide <- dcast.data.table(sc2.stm[aeid==1691 & !dsstox_substance_id %in% stm.mc.dtxsids & hitc==0], dsstox_substance_id + chnm ~aenm,
                             value.var='hitc',
                             fun.aggregate=max)

add.wide[,STM_H9_Viability_norm := NA]

#length(unique(add.wide$dsstox_substance_id)) #133

3.3 Calculate dev flag

• Calculate the DEV flag.
• Any hit in STM_H9_OrnCyssISnorm_ratio_dn considered relevant as a DEV flag.
• Hits where the distance from the parallel viability assay is more then 0.25 log are considered selective flags.

dev.flag <- rbind(mc5.stm.wide,
                  add.wide, fill=TRUE)

dev.flag[,stm.cyto.dist := ifelse(!is.na(STM_H9_Viability_norm), STM_H9_Viability_norm - STM_H9_OrnCyssISnorm_ratio_dn, 3 - STM_H9_OrnCyssISnorm_ratio_dn)]

dev.flag[,dev.flag := 0]
dev.flag[!is.na(STM_H9_OrnCyssISnorm_ratio_dn), dev.flag := 1]
dev.flag[,dev.flag.specific := 0]
dev.flag[stm.cyto.dist >0.25, dev.flag.specific := 1]

Add TEST DEV TOX prediction and the applicability domain.

load('./source/chem/apcra_chem_ad.RData')
setnames(test.opera.pred.total, c('DTXSID','PREFERRED_NAME'), c('dsstox_substance_id','chnm'), skip_absent = TRUE)
dev.flag$test.devtox.score <- test.opera.pred.total$DEVTOX_TEST_PRED[match(dev.flag$dsstox_substance_id,
                                                               test.opera.pred.total$dsstox_substance_id)]

# replace Inf with NA
for (j in 1:ncol(dev.flag)) set(dev.flag, which(is.infinite(dev.flag[[j]])), j, NA)
# replace - with NA
dev.flag[,c('test.devtox.score')] <- lapply(dev.flag[,c('test.devtox.score')], function(col) as.numeric(gsub("-$|\\,",NA, col)))
dev.flag$aqc_indicator <- ad.tbl$aqc_indicator[match(dev.flag$dsstox_substance_id, ad.tbl$DTXSID)]
dev.flag$apcra.pro.only <- ad.tbl$apcra.pro.only[match(dev.flag$dsstox_substance_id, ad.tbl$DTXSID)]

3.4 Create Endocrine Hazard Flags

• Match the logic in the recent comments provided to the EPA program office partners on a prioritization and screening workflow using ER and AR models.
  • If ToxCast ER/AR pathway models available: these scores trump CERAPP and COMPARA.
  • ToxCast ER/AR pathway model positives: => 0.1
  • Grouped equivocals (0.001-0.1) with the negative to enable easier prioritization.
  • For ToxCast AR pathway model (older version to match what is in invitrodb version 3.3 and CompTox Chemicals Dashboard), require confidence flags to be > 2 for a positive.
  • Highlight substances with no data - could be a data gap?

total.endo.dtxsid <- c(apcra.list$DTXSID, dart.ref)

## here we should really get the dtxsids for cerapp...

cerapp.compara <- merge.data.table(cerapp[,c('CASRN','CHEMICAL_NAME','input_SMILES','InChI_Key',
                                                'Potency_class_2_binding', 'Potency_class_2_agonist', 'Potency_class_2_antagonist','consensus_2_binding', 'consensus_2_agonist','consensus_2_antagonist')],
                                by.x = 'CASRN',
                                all.x=TRUE,
                                compara[, c('dsstox_substance_id','casrn','consensus_binding', 'consensus_agonist','consensus_antagonist')],
                                by.y='casrn',
                                all.y=TRUE)

colnames(er) <- paste0('er_model_', colnames(er))
er.cerapp.compara <- merge.data.table(cerapp.compara,
                                er,
                                by.x='CASRN',
                                by.y='er_model_CASRN',
                                all.x=TRUE)

colnames(ar) <- paste0('ar_model_', colnames(ar))
er.ar.cerapp.compara <- merge.data.table(er.cerapp.compara,
                                ar,
                                by.x='CASRN',
                                by.y='ar_model_CASRN',
                                all.x=TRUE)

er.ar.apcra <- er.ar.cerapp.compara[dsstox_substance_id %in% total.endo.dtxsid]

col.num <- c('consensus_2_binding', 'consensus_2_agonist', 'consensus_2_antagonist')
er.ar.apcra[, (col.num) := lapply (.SD, as.numeric), .SDcols = col.num ]
# indicate any active result from compara consensus qsar
er.ar.apcra[is.na(consensus_binding & consensus_agonist & consensus_antagonist),compara.flag := 2] # if data not available
er.ar.apcra[consensus_binding==0 & consensus_agonist==0 & consensus_antagonist==0, compara.flag := 0] # if data are all negative
er.ar.apcra[consensus_binding==1|consensus_agonist==1|consensus_antagonist==1, compara.flag := 1] # if any mode is positive

# indicate any active result from cerapp consensus qsar
er.ar.apcra[is.na(consensus_2_agonist & consensus_2_binding & consensus_2_antagonist),cerapp.flag := 2]
er.ar.apcra[consensus_2_agonist==0 & consensus_2_binding==0 & consensus_2_antagonist==0, cerapp.flag := 0]
er.ar.apcra[consensus_2_agonist==1|consensus_2_binding==1|consensus_2_antagonist==1, cerapp.flag := 1]

# create positive, equivocal, and negative code on ToxCast AR model
er.ar.apcra[is.na(ar_model_Agonist_AUC & ar_model_Antagonist_AUC), toxcast.ar.flag := 2] # data not available
er.ar.apcra[ar_model_Agonist_AUC < 0.1 & ar_model_Antagonist_AUC < 0.1|ar_model_Antagonist_AUC > 0.1 & ar_model_Antagonist_Confidence_Score <= 2, toxcast.ar.flag := 0] # grouped equivocals into the negative space
er.ar.apcra[ar_model_Agonist_AUC >= 0.1|ar_model_Antagonist_AUC >= 0.1 & ar_model_Antagonist_Confidence_Score >2, 
            toxcast.ar.flag := 1]

# create positive, equivocal, and negative code on ToxCast ER model
er.ar.apcra[is.na(er_model_Agonist_AUC & er_model_Antagonist_AUC), toxcast.er.flag := 2] # no data available
er.ar.apcra[er_model_Agonist_AUC< 0.1 & er_model_Antagonist_AUC < 0.1, toxcast.er.flag := 0] # grouped equivocals into the negative space
er.ar.apcra[er_model_Agonist_AUC >= 0.1|er_model_Antagonist_AUC >= 0.1, toxcast.er.flag := 1]

# require toxcast er model data to trump cerapp call
er.ar.apcra[toxcast.er.flag ==2 & cerapp.flag==2, final.er.flag :=0.1] # no data
er.ar.apcra[toxcast.er.flag ==2 & cerapp.flag==0, final.er.flag :=0] # only cerapp data and it's negative
er.ar.apcra[toxcast.er.flag ==2 & cerapp.flag==1, final.er.flag :=0.5] # only cerapp data and it's positive
er.ar.apcra[toxcast.er.flag ==0 & cerapp.flag %in% c(0,1,2), final.er.flag := 0] # er model available and negative
er.ar.apcra[toxcast.er.flag ==1 & cerapp.flag %in% c(0,1,2), final.er.flag := 1] # er model available and positive

# require toxcast ar model data to trump compara call
er.ar.apcra[toxcast.ar.flag ==2 & compara.flag==2, final.ar.flag :=0.1] # no data
er.ar.apcra[toxcast.ar.flag ==2 & compara.flag==0, final.ar.flag :=0] # only cerapp data and it's negative
er.ar.apcra[toxcast.ar.flag ==2 & compara.flag==1, final.ar.flag :=0.5] # only cerapp data and it's positive
er.ar.apcra[toxcast.ar.flag ==0 & compara.flag %in% c(0,1,2), final.ar.flag := 0] # er model available and negative
er.ar.apcra[toxcast.ar.flag ==1 & compara.flag %in% c(0,1,2), final.ar.flag := 1] # er model available and positive


er.ar.apcra[, ed.flag.sum := sum(final.er.flag, final.ar.flag), by=list(dsstox_substance_id)]

er.ar.flag <- er.ar.apcra[,c("dsstox_substance_id","compara.flag","cerapp.flag","toxcast.ar.flag","toxcast.er.flag","final.er.flag",                          "final.ar.flag","ed.flag.sum")]

3.5 DART flag

load('./output/APCRA_haz_flg_BER.RData')

# add minimum NAM column for comparison
dev.flag$min_nam <- asnm.tier1$col_min[match(dev.flag$dsstox_substance_id,
                                     asnm.tier1$dsstox_substance_id)]
dev.flag[min_nam=='CCTE', min_nam := 'MEA'] # the CCTE vendor is for the acute MEA data from Shafer lab

# define the binary flag for the TEST dev model
dev.flag[test.devtox.score > 0.7, dev.test := 0.5]
dev.flag[test.devtox.score < 0.7, dev.test := 0]

# merge dev and ER/AR flag
dart.flag <- merge.data.table(dev.flag,
                              er.ar.flag,
                              by=c('dsstox_substance_id'),
                              all.x=TRUE,
                              all.y=TRUE)

# make sure chnm is not missing
setnames(apcra.total, 'DTXSID', 'dsstox_substance_id')
setnames(apcra.total, 'preferred_name', 'chnm')
dart.flag.na <- dart.flag[is.na(chnm)]
dart.flag$chnm <- apcra.total$chnm[match(dart.flag$dsstox_substance_id,
                                            apcra.total$dsstox_substance_id)]
dart.flag$chnm_dev <- dev.flag$chnm[match(dart.flag$dsstox_substance_id,
                                          dev.flag$dsstox_substance_id)]
dart.flag[is.na(chnm), chnm := chnm_dev]
dart.flag[,c('chnm_dev') := NULL]

# add BER
dart.flag$ber.med.aed50 <- asnm.tier1.all.invivo$ber.med.aed50[match(dart.flag$dsstox_substance_id,
                                                                     asnm.tier1.all.invivo$dsstox_substance_id)]

dart.flag <- dart.flag[order(ber.med.aed50)]

# further refine columns of data available
dart_mat <- dart.flag[apcra.pro.only==1]
dart_mat <- dart_mat[aqc_indicator==1]
dart_mat <- dart_mat[ber.med.aed50 < 4 ]
dart_ref_mat <- unique(dart.flag[dsstox_substance_id %in% dart.ref])
dart_ref_mat$chnm <- dev.flag$chnm[match(dart_ref_mat$dsstox_substance_id,
                                         dev.flag$dsstox_substance_id)]

setnames(dart_mat, c('dev.test',
                                       'dev.flag',
                                       'dev.flag.specific',
                                       'final.er.flag',
                                       'final.ar.flag'),
         c('DEV-TEST','DEV','DEV-S','ER','AR'))
setnames(dart_ref_mat, c('dev.test',
                                       'dev.flag',
                                       'dev.flag.specific',
                                       'final.er.flag',
                                       'final.ar.flag'),
         c('DEV-TEST','DEV','DEV-S','ER','AR'))

# comprise the main matrices
dart_mat2 <- as.matrix(dart_mat[,c('DEV-TEST','DEV','DEV-S','ER','AR')])
rownames(dart_mat2) <- dart_mat[, chnm] # define rownames as chemical name

dart_ref2 <- as.matrix(dart_ref_mat[,c('DEV-TEST','DEV','DEV-S','ER','AR')])
rownames(dart_ref2) <- dart_ref_mat[, chnm] # define rownames as chemical name

head(dart_ref2)

##                    DEV-TEST DEV DEV-S ER AR
## Flutamide               0.5   1     1  0  1
## Boric acid (H3BO3)       NA   0     0 NA NA
## 5-Fluorouracil           NA   1     0  0  0
## Diethylstilbestrol      0.5   0     0  1  1
## Vinclozolin             0.5   0     0  0  1
## Genistein               0.5   1     0  1  0

# annotations for main hm
anno_df <- data.frame(rownames(dart_mat2))
anno_df$min_nam <- dart.flag$min_nam[match(anno_df[,1],
                                            dart.flag$chnm)]
anno_df$chnm <- dart.flag$chnm[match(anno_df[,1],
                                      dart.flag$chnm)]
anno_df$ber <- dart.flag$ber.med.aed50[match(anno_df[,1],
                                             dart.flag$chnm)]

#left_ha <- rowAnnotation(min_nam = anno_text(anno_df$min_nam,
#                                            just='right',
#                                            location=1,
#                                            show_name = FALSE))
right_ha <- rowAnnotation(BER = anno_barplot(anno_df$ber, width = unit(4,'cm')))

# annotations for ref hm
anno_ref <- data.frame(rownames(dart_ref2))
anno_ref$chnm <- dart.flag$chnm[match(anno_ref[,1],
                                      dart.flag$chnm)]

hm_dart <- Heatmap(matrix = dart_mat2, 
                       cluster_columns = FALSE,
                       cluster_rows=FALSE,
                       name="log10-mg/kg/day",
                       #col=col_fun,
                       na_col='gray',
                       #col= colors,
                       col = colorRamp2(breaks = c(0, 0.5, 1),
                                         colors = c("white" , "#3CBB75FF", "#440154FF")),
                       #col = list(type=c("1" = "#2166ac","0"= "#f7f7f7")),
                       show_row_names = TRUE, 
                       row_dend_width = unit(3, "cm"),
                       column_names_max_height = unit(8, "cm"),
                       column_names_gp = gpar(fontsize = 12),
                       heatmap_legend_param = list(title = "In silico/In vitro",legend_direction = 'horizontal'),
                       width=unit(6,'cm'),
                       height = unit(14, 'cm'),
                       rect_gp = gpar(col = "gray", lwd = 2),
                       column_names_side='top',
                       #left_annotation = left_ha,
                       right_annotation = right_ha,
                   row_names_gp=gpar(fontsize=10,fontfamily='sans'))

hm_dart_ref <- Heatmap(matrix = dart_ref2, 
                       cluster_columns = FALSE,
                       cluster_rows=FALSE,
                       name="log10-mg/kg/day",
                       #col=col_fun,
                       na_col='gray',
                       #col= colors,
                       col = colorRamp2(breaks = c(0, 0.5, 1),
                                         colors = c("white" , "#3CBB75FF", "#440154FF")),
                       #col = list(type=c("1" = "#2166ac","0"= "#f7f7f7")),
                       show_row_names = TRUE, 
                       #row_dend_width = unit(3, "cm"),
                       #column_names_max_height = unit(8, "cm"),
                       column_names_gp = gpar(fontsize = 12),
                       heatmap_legend_param = list(title = "In silico/In vitro",legend_direction = 'horizontal'),
                       #width=unit(6,'cm'),
                       height = unit(4, 'cm'),
                       rect_gp = gpar(col = "gray", lwd = 2),
                       column_names_side='top',
                   row_names_gp=gpar(fontsize=10,fontfamily='sans'))

hm.ed <- draw(hm_dart %v%  hm_dart_ref, 
              heatmap_legend_side='bottom',
              ht_gap = unit(1, "cm"))

file.dir <- paste("./output", sep="")
file.name <- paste("/Figure_DART_HM_", Sys.Date(), ".png", sep="")
file.path <- paste(file.dir, file.name, sep="")
dir.create(path=file.dir, showWarnings = FALSE, recursive = TRUE)
png(file.path, width = 4000, height = 5000, res=450)
hm.ed
dev.off()

## png 
##   2

4 BSK data

The BSK data, or BioMAP panel as it is now known, contain endpoints potentially relevant to immunosuppression. We need to annotate the specific endpoints to calculate the potency of potential immunosuppression and the potency of potential “acute toxicity” (or really, cytotoxicity) in these pathophysiological systems.

4.1 Annotate the BSK assay endpoints

bsk.aeid <- unique(mega.mc5[asnm=='BSK',c('aeid','aenm')])
bsk.aeid[,cells := tstrsplit(aenm, "_", fixed=TRUE, keep=c(2))]
bsk.aeid[,endpoint := tstrsplit(aenm, "_", fixed=TRUE, keep=c(3))]
bsk.aeid[,direction := tstrsplit(aenm, "_", fixed=TRUE, keep=c(4))]

## label signature activity groups
bsk.aeid[,activity := as.character(NA)]
bsk.aeid[grep('SRB_down', aenm),activity := 'cytotoxicity']
bsk.aeid[grep('_3C_Proliferation_down',aenm), activity := 'liver toxicity'] #endothelial cell signal
bsk.aeid[grep('SAg_Proliferation_down',aenm), activity := 'immunosuppression'] # decreased T-cells
bsk.aeid[grep('BT_sIgG_down', aenm), activity := 'immunosuppression'] # decreased IgG
bsk.aeid[grep('BT_Bcell_Proliferation_down',aenm), activity := 'immunosuppression'] # decreased B cell proliferation
bsk.aeid[grep('PBMCCytotoxicity_down', aenm), activity := 'immunosuppression'] #cytotoxic to peripheral blood mononuclear cells
bsk.aeid[grep('3C_Thrombomodulin_up', aenm), activity := 'thrombosis'] 

# next systems are at 'non-cytotoxic concentrations'
bsk.aeid[grep('LPS_PGE2_up', aenm), activity := 'skin irritation']
bsk.aeid[grep('LPS_TNFa_up', aenm), activity := 'skin irritation']
bsk.aeid[grep('hDFCGF_CollagenIII_down', aenm), activity := 'skin sensitization']
bsk.aeid[grep('hDFCGF_VCAM1_up', aenm), activity := 'skin rash']
bsk.aeid[grep('CASM3C_SAA_up', aenm), activity := 'vascular toxicity']
bsk.aeid[aenm=='BSK_IMphg_IL10_down', activity := 'immunosuppression']

kable(bsk.aeid[activity %in% c('immunosuppression','skin irritation','skin rash','skin sensitization')])

aeid	aenm	cells	endpoint	direction	activity
235	BSK_hDFCGF_CollagenIII_down	hDFCGF	CollagenIII	down	skin sensitization
258	BSK_hDFCGF_VCAM1_up	hDFCGF	VCAM1	up	skin rash
290	BSK_LPS_PGE2_up	LPS	PGE2	up	skin irritation
296	BSK_LPS_TNFa_up	LPS	TNFa	up	skin irritation
313	BSK_SAg_PBMCCytotoxicity_down	SAg	PBMCCytotoxicity	down	immunosuppression
315	BSK_SAg_Proliferation_down	SAg	Proliferation	down	immunosuppression
2810	BSK_BT_Bcell_Proliferation_down	BT	Bcell	Proliferation	immunosuppression
2812	BSK_BT_PBMCCytotoxicity_down	BT	PBMCCytotoxicity	down	immunosuppression
2814	BSK_BT_sIgG_down	BT	sIgG	down	immunosuppression
2928	BSK_IMphg_IL10_down	IMphg	IL10	down	immunosuppression

4.2 Load positive controls for immunosuppression

In the BioMAP (BSK) panel, we also ran immunosuppressive drugs as controls. We can use these as reference chemicals, but they are not in the APCRA list of chemicals and so we need to retrieve this from invitrodb v3.5.

Note that in invitrodb v3.5, BSK data are represented as ac50 (modl_ga), but are actually lowest effect concentrations. This modeling choice was related to the low number of concentrations (most often 4), low number of replicates (most often 2), and low top-over-cutoff typically observed.Standard curve-fitting models had resulted in extremely low ac50s and so we chose a lowest effect concentration approach.

immunosupp.drug.dtxsids <- c('DTXSID0020365', # cyclosporin A
                             'DTXSID3047429', # dexamethasone sodium phosphate
                             'DTXSID4020119', # azathioprine
                             'DTXSID4020822' # methotrexate
                             ) 

spids <- tcplLoadChem(field='dsstox_substance_id',val=immunosupp.drug.dtxsids) # get sample ids
mc5.bsk.ref <- tcplPrepOtpt(tcplLoadData(lvl=5,type='mc',fld='aeid', val=bsk.aeid$aeid)) # get mc5 data
mc5.bsk.ref <- mc5.bsk.ref[dsstox_substance_id %in% spids$dsstox_substance_id] # narrow mc5 data to immunosuppressive drugs

mc5.bsk.combined <- rbind(mega.mc5[asnm=='BSK'],
                          mc5.bsk.ref,
                          fill=TRUE)

bsk.flag <- unique(mc5.bsk.combined[ , list(
  total.endpoints.screened  = .N, #total number of aeids tested in mc
  active.assay.count  = as.double(length(which(hitc==1))),  # active count
  inactive.assay.count  = as.double(length(which(hitc==0))),  #inactive count
  active.percent = round((length(which(hitc==1))/.N)*100,2), #active percent
  inactive.percent = round((length(which(hitc==0))/.N)*100,2), #inactive percent
  min.log.lec = min(modl_ga, na.rm=TRUE), 
  med.log.lec = median(modl_ga, na.rm=TRUE),
  acute.tox = min(modl_ga[aenm %in% c("BSK_3C_SRB_down", 
                                      "BSK_4H_SRB_down", 
                                      "BSK_BE3C_SRB_down", 
                                      "BSK_CASM3C_SRB_down", 
                                      "BSK_hDFCGF_SRB_down",
                                      "BSK_KF3CT_SRB_down", 
                                      "BSK_LPS_SRB_down", 
                                      "BSK_SAg_SRB_down", 
                                      "BSK_MyoF_SRB_down", 
                                      "BSK_BF4T_SRB_down", 
                                      "BSK_IMphg_SRB_down")], na.rm=TRUE),
  skin.irrit.lec = min(modl_ga[aenm=='BSK_LPS_PGE2_up'], na.rm=TRUE),
  skin.irrit.spec = (min(modl_ga[aenm %in% c("BSK_3C_SRB_down", 
                                             "BSK_4H_SRB_down", 
                                             "BSK_BE3C_SRB_down", 
                                             "BSK_CASM3C_SRB_down", 
                                             "BSK_hDFCGF_SRB_down", 
                                             "BSK_KF3CT_SRB_down", 
                                             "BSK_LPS_SRB_down", 
                                             "BSK_SAg_SRB_down", 
                                             "BSK_MyoF_SRB_down", 
                                             "BSK_BF4T_SRB_down", 
                                             "BSK_IMphg_SRB_down")], na.rm=TRUE)) - (min(modl_ga[aenm=='BSK_LPS_PGE2_up'], na.rm=TRUE)), 
  skin.rash.lec = min(modl_ga[aenm %in% c('BSK_hDFCGF_VCAM1_up')], na.rm=TRUE),
  skin.rash.spec = (min(modl_ga[aenm %in% c("BSK_3C_SRB_down", "BSK_4H_SRB_down", "BSK_BE3C_SRB_down", 
"BSK_CASM3C_SRB_down", "BSK_hDFCGF_SRB_down", "BSK_KF3CT_SRB_down", 
"BSK_LPS_SRB_down", "BSK_SAg_SRB_down", "BSK_MyoF_SRB_down", 
"BSK_BF4T_SRB_down", "BSK_IMphg_SRB_down")], na.rm=TRUE)) - (min(modl_ga[aenm %in% c('BSK_hDFCGF_VCAM1_up')], na.rm=TRUE)),
  skin.sens.lec = min(modl_ga[aenm=='BSK_hDFCGF_CollagenIII_down']),
  skin.sens.spec = (min(modl_ga[aenm %in% c("BSK_3C_SRB_down", "BSK_4H_SRB_down", "BSK_BE3C_SRB_down", 
"BSK_CASM3C_SRB_down", "BSK_hDFCGF_SRB_down", "BSK_KF3CT_SRB_down", 
"BSK_LPS_SRB_down", "BSK_SAg_SRB_down", "BSK_MyoF_SRB_down", 
"BSK_BF4T_SRB_down", "BSK_IMphg_SRB_down")], na.rm=TRUE)) - (min(modl_ga[aenm %in% c('BSK_hDFCGF_CollagenIII_down')], na.rm=TRUE)),
skin.inflamm.lec = min(modl_ga[aenm %in% c('BSK_hDFCGF_CollagenIII_down','BSK_hDFCGF_VCAM1_up','BSK_LPS_PGE2_up')], na.rm=TRUE),
skin.inflamm.spec = (min(modl_ga[aenm %in% c("BSK_3C_SRB_down", "BSK_4H_SRB_down", "BSK_BE3C_SRB_down", 
"BSK_CASM3C_SRB_down", "BSK_hDFCGF_SRB_down", "BSK_KF3CT_SRB_down", 
"BSK_LPS_SRB_down", "BSK_SAg_SRB_down", "BSK_MyoF_SRB_down", 
"BSK_BF4T_SRB_down", "BSK_IMphg_SRB_down")], na.rm=TRUE)) - (min(modl_ga[aenm %in% c('BSK_hDFCGF_CollagenIII_down','BSK_hDFCGF_VCAM1_up','BSK_LPS_PGE2_up')], na.rm=TRUE)),
  immun.sup.lec = min(modl_ga[aenm %in% bsk.aeid[activity=='immunosuppression']$aenm], na.rm=TRUE),
  immun.sup.spec = (min(modl_ga[aenm %in% c("BSK_3C_SRB_down", "BSK_4H_SRB_down", "BSK_BE3C_SRB_down", 
"BSK_CASM3C_SRB_down", "BSK_hDFCGF_SRB_down", "BSK_KF3CT_SRB_down", 
"BSK_LPS_SRB_down", "BSK_SAg_SRB_down", "BSK_MyoF_SRB_down", 
"BSK_BF4T_SRB_down", "BSK_IMphg_SRB_down")], na.rm=TRUE)) - (min(modl_ga[aenm %in% c('BSK_3C_Proliferation_down')], na.rm=TRUE))
), by = list(dsstox_substance_id, chnm, casn)]) # could make this by spid because if there are n>1 spid the assay number will mess up

invisible(lapply(names(bsk.flag), function(.name) set(bsk.flag, which(is.infinite(bsk.flag[[.name]])), j=.name, value=NA))) # remove Inf
invisible(lapply(names(bsk.flag), function(.name) set(bsk.flag, which(is.nan(bsk.flag[[.name]])), j=.name, value=NA))) # remove NaN

bsk.flag[,selective := ifelse(!is.na(acute.tox), acute.tox-immun.sup.lec, 3-immun.sup.lec)]
colnames(bsk.flag)

Wrangling the data to visualize and adding the applicability domain information.

bsk.flag2 <- bsk.flag[,c('dsstox_substance_id',
                        'chnm',
                        'acute.tox',
                        'immun.sup.lec',
                        'selective')]

bsk.flag2$aqc_indicator <- ad.tbl$aqc_indicator[match(bsk.flag2$dsstox_substance_id, ad.tbl$DTXSID)]
bsk.flag2$apcra.pro.only <- ad.tbl$apcra.pro.only[match(bsk.flag2$dsstox_substance_id, ad.tbl$DTXSID)]

# add BER
bsk.flag2$ber.med.aed50 <- asnm.tier1.all.invivo$ber.med.aed50[match(bsk.flag2$dsstox_substance_id,
                                                                     asnm.tier1.all.invivo$dsstox_substance_id)]

bsk.flag2 <- bsk.flag2[order(ber.med.aed50)]

5 MEA data

Examine the pattern of the MEA data. In the end this was not very fruitful.

5.1 Define activity type for the aeid’s in the MEA acute.

mea.aeid <- unique(mega.mc5[asnm=='CCTE', c('aeid','aenm')])
mea.aeid[,process := tstrsplit(aenm, "CCTE_Shafer_MEA_", keep=c(2))]
mea.aeid[, direction := str_extract(aenm, "[a-z]{2}$")]

## label signature activity groups

### General Firing Activity
mea.aeid[,activity := as.character(NA)]
mea.aeid[grep('firing', aenm),activity := 'Firing']
mea.aeid[grep('MFR',aenm), activity := 'Firing']
mea.aeid[grep('acute_burst_number',aenm), activity := "Firing"]
mea.aeid[grep('acute_spike_number',aenm), activity := "Firing"]

### Burst structure
mea.aeid[is.na(activity) & grep('burst',aenm), activity := "Bursting"]
mea.aeid[is.na(activity) & grep('spike',aenm), activity := "Bursting"]

### Connectivity
mea.aeid[grep('network', aenm), activity := 'Connectivity']
mea.aeid[grep('cross_correlation', aenm), activity := 'Connectivity']
mea.aeid[grep('synchrony', aenm), activity := 'Connectivity']

### Cytotoxicity
mea.aeid[grep('LDH', aenm), activity := 'Cytotoxicity']
mea.aeid[grep('AB', aenm), activity := 'Cytotoxicity']

mea.aeid <- mea.aeid[order(activity,aeid)]
# write.csv
#write.csv(mea.aeid,"./output/mea_aeid_annotation_17apr2021.csv")

# examine
kable(mea.aeid, 
          filter='top', 
          options=list(pagelength=25, autoWidth=FALSE,  scrollX=TRUE, initComplete = JS(
    "function(settings, json) {",
    "$('body').css({'font-family': 'Calibri'});",
    "}"
  )))

aeid	aenm	process	direction	activity
2460	CCTE_Shafer_MEA_acute_burst_duration_mean_dn	acute_burst_duration_mean_dn	dn	Bursting
2461	CCTE_Shafer_MEA_acute_burst_duration_mean_up	acute_burst_duration_mean_up	up	Bursting
2462	CCTE_Shafer_MEA_acute_per_burst_spike_number_mean_dn	acute_per_burst_spike_number_mean_dn	dn	Bursting
2463	CCTE_Shafer_MEA_acute_per_burst_spike_number_mean_up	acute_per_burst_spike_number_mean_up	up	Bursting
2464	CCTE_Shafer_MEA_acute_interburst_interval_mean_dn	acute_interburst_interval_mean_dn	dn	Bursting
2465	CCTE_Shafer_MEA_acute_interburst_interval_mean_up	acute_interburst_interval_mean_up	up	Bursting
2466	CCTE_Shafer_MEA_acute_burst_percentage_mean_dn	acute_burst_percentage_mean_dn	dn	Bursting
2467	CCTE_Shafer_MEA_acute_burst_percentage_mean_up	acute_burst_percentage_mean_up	up	Bursting
2468	CCTE_Shafer_MEA_acute_burst_percentage_std_dn	acute_burst_percentage_std_dn	dn	Bursting
2469	CCTE_Shafer_MEA_acute_burst_percentage_std_up	acute_burst_percentage_std_up	up	Bursting
2470	CCTE_Shafer_MEA_acute_per_network_burst_spike_number_mean_dn	acute_per_network_burst_spike_number_mean_dn	dn	Connectivity
2471	CCTE_Shafer_MEA_acute_per_network_burst_spike_number_mean_up	acute_per_network_burst_spike_number_mean_up	up	Connectivity
2472	CCTE_Shafer_MEA_acute_per_network_burst_spike_number_std_dn	acute_per_network_burst_spike_number_std_dn	dn	Connectivity
2473	CCTE_Shafer_MEA_acute_per_network_burst_spike_number_std_up	acute_per_network_burst_spike_number_std_up	up	Connectivity
2474	CCTE_Shafer_MEA_acute_per_network_burst_electrodes_number_mean_dn	acute_per_network_burst_electrodes_number_mean_dn	dn	Connectivity
2475	CCTE_Shafer_MEA_acute_per_network_burst_electrodes_number_mean_up	acute_per_network_burst_electrodes_number_mean_up	up	Connectivity
2476	CCTE_Shafer_MEA_acute_network_burst_percentage_dn	acute_network_burst_percentage_dn	dn	Connectivity
2477	CCTE_Shafer_MEA_acute_network_burst_percentage_up	acute_network_burst_percentage_up	up	Connectivity
2478	CCTE_Shafer_MEA_acute_cross_correlation_area_dn	acute_cross_correlation_area_dn	dn	Connectivity
2479	CCTE_Shafer_MEA_acute_cross_correlation_area_up	acute_cross_correlation_area_up	up	Connectivity
2480	CCTE_Shafer_MEA_acute_cross_correlation_HWHM_dn	acute_cross_correlation_HWHM_dn	dn	Connectivity
2481	CCTE_Shafer_MEA_acute_cross_correlation_HWHM_up	acute_cross_correlation_HWHM_up	up	Connectivity
2482	CCTE_Shafer_MEA_acute_synchrony_index_dn	acute_synchrony_index_dn	dn	Connectivity
2483	CCTE_Shafer_MEA_acute_synchrony_index_up	acute_synchrony_index_up	up	Connectivity
2540	CCTE_Shafer_MEA_acute_LDH_up	acute_LDH_up	up	Cytotoxicity
2541	CCTE_Shafer_MEA_acute_AB_dn	acute_AB_dn	dn	Cytotoxicity
2454	CCTE_Shafer_MEA_acute_spike_number_dn	acute_spike_number_dn	dn	Firing
2455	CCTE_Shafer_MEA_acute_spike_number_up	acute_spike_number_up	up	Firing
2456	CCTE_Shafer_MEA_acute_firing_rate_mean_dn	acute_firing_rate_mean_dn	dn	Firing
2457	CCTE_Shafer_MEA_acute_firing_rate_mean_up	acute_firing_rate_mean_up	up	Firing
2458	CCTE_Shafer_MEA_acute_burst_number_dn	acute_burst_number_dn	dn	Firing
2459	CCTE_Shafer_MEA_acute_burst_number_up	acute_burst_number_up	up	Firing

5.2 Make MEA flag

mea.flag <- unique(mega.mc5[asnm=='CCTE' , list(
  total.endpoints.screened  = .N, #total number of aeids tested in mc
  active.assay.count  = as.double(length(which(hitc==1))),  # active count
  inactive.assay.count  = as.double(length(which(hitc==0))),  #inactive count
  active.percent = round((length(which(hitc==1))/.N)*100,2), #active percent
  inactive.percent = round((length(which(hitc==0))/.N)*100,2), #inactive percent
  min.log.acc = min(modl_acc, na.rm=TRUE), 
  med.log.acc = median(modl_acc, na.rm=TRUE),
  p5.log.acc = quantile(modl_acc, c(0.05),na.rm=TRUE),
  up.ct = as.double(length(which(hitc==1 & aeid %in% mea.aeid[direction=='up']$aeid))),
  dn.ct = as.double(length(which(hitc==1 & aeid %in% mea.aeid[direction=='dn']$aeid))),
  bursting = min(modl_acc[aeid %in% mea.aeid[activity=='Bursting']$aeid], na.rm=TRUE),
  connectivity = min(modl_acc[aeid %in% mea.aeid[activity=='Connectivity']$aeid], na.rm=TRUE),
  firing = min(modl_acc[aeid %in% mea.aeid[activity=='Firing']$aeid], na.rm=TRUE),
  cytotoxicity = min(modl_acc[aeid %in% mea.aeid[activity=='Cytotoxicity']$aeid], na.rm=TRUE)
  ), by = list(dsstox_substance_id, chnm, casn, spid)]) # could make this by spid because if there are n>1 spid the assay number will mess up

invisible(lapply(names(mea.flag), function(.name) set(mea.flag, which(is.infinite(mea.flag[[.name]])), j=.name, value=NA))) # remove Inf
invisible(lapply(names(mea.flag), function(.name) set(mea.flag, which(is.nan(mea.flag[[.name]])), j=.name, value=NA))) # remove NaN

mea.flag[,bursting.spec := ifelse(!is.na(cytotoxicity), cytotoxicity - bursting, 3-bursting), by=list(dsstox_substance_id)]
mea.flag[,connectivity.spec := ifelse(!is.na(cytotoxicity), cytotoxicity - connectivity, 3-connectivity), by=list(dsstox_substance_id)]
mea.flag[,firing.spec := ifelse(!is.na(cytotoxicity), cytotoxicity - firing, 3-firing), by=list(dsstox_substance_id)]

• get some reference chemicals

mea <- tcplPrepOtpt(tcplLoadData(lvl=5,type='mc',fld='aeid',val=mea.aeid$aeid))
#mea.ref <- mea[ chnm %in% c('Tributyltin chloride','Abamectin','Lindane','beta-Cyfluthrin')]

mea.ref <- mea[dsstox_substance_id %in% c(
  'DTXSID2020686', #lindane
  'DTXSID3027403', # tributyltin chloride
  'DTXSID8023892', # abamectin
  'DTXSID8032330', # beta-Cyfluthrin
  'DTXSID9058238')] # avermectin B1a

mea.ref.dtxsids <- c(
  'DTXSID2020686', #lindane
  'DTXSID3027403', # tributyltin chloride
  'DTXSID8023892', # abamectin
  'DTXSID8032330', # beta-Cyfluthrin
  'DTXSID9058238') # avermectin B1a

mea.flag.ref <- unique(mea.ref[ , list(
  total.endpoints.screened  = .N, #total number of aeids tested in mc
  active.assay.count  = as.double(length(which(hitc==1))),  # active count
  inactive.assay.count  = as.double(length(which(hitc==0))),  #inactive count
  active.percent = round((length(which(hitc==1))/.N)*100,2), #active percent
  inactive.percent = round((length(which(hitc==0))/.N)*100,2), #inactive percent
  min.log.acc = min(modl_acc, na.rm=TRUE), 
  med.log.acc = median(modl_acc, na.rm=TRUE),
  p5.log.acc = quantile(modl_acc, c(0.05),na.rm=TRUE),
  up.ct = as.double(length(which(hitc==1 & aeid %in% mea.aeid[direction=='up']$aeid))),
  dn.ct = as.double(length(which(hitc==1 & aeid %in% mea.aeid[direction=='dn']$aeid))),
  bursting = min(modl_acc[aeid %in% mea.aeid[activity=='Bursting']$aeid], na.rm=TRUE),
  connectivity = min(modl_acc[aeid %in% mea.aeid[activity=='Connectivity']$aeid], na.rm=TRUE),
  firing = min(modl_acc[aeid %in% mea.aeid[activity=='Firing']$aeid], na.rm=TRUE),
  cytotoxicity = min(modl_acc[aeid %in% mea.aeid[activity=='Cytotoxicity']$aeid], na.rm=TRUE)
  ), by = list(dsstox_substance_id, chnm, casn)]) # could make this by spid because if there are n>1 spid the assay number will mess up

invisible(lapply(names(mea.flag.ref), function(.name) set(mea.flag.ref, which(is.infinite(mea.flag.ref[[.name]])), j=.name, value=NA))) # remove Inf
invisible(lapply(names(mea.flag.ref), function(.name) set(mea.flag.ref, which(is.nan(mea.ref[[.name]])), j=.name, value=NA))) # remove NaN

mea.flag.ref[,bursting.spec := ifelse(!is.na(cytotoxicity), cytotoxicity - bursting, 3-bursting), by=list(dsstox_substance_id)]
mea.flag.ref[,connectivity.spec := ifelse(!is.na(cytotoxicity), cytotoxicity - connectivity, 3-connectivity), by=list(dsstox_substance_id)]
mea.flag.ref[,firing.spec := ifelse(!is.na(cytotoxicity), cytotoxicity - firing, 3-firing), by=list(dsstox_substance_id)]

Importantly, we only let the MEA flag be active if > 3 endpoints in the same direction are positive. The 5th percentile ACC of the positive MEA chemicals will serve as a potency comparator/flag (if most sensitive potency for the chemical).

mea.flag2 <- rbind(mea.flag[,-c('spid')], mea.flag.ref, fill=TRUE)

#mea.flag2[dsstox_substance_id %in% c(
#  'DTXSID2020686', #lindane
#  'DTXSID3027403', # tributyltin chloride - clearly tested in 2 spids
#  'DTXSID8023892', # abamectin
#  'DTXSID8032330', # beta-Cyfluthrin
#  'DTXSID9058238')]

mea.flag3 <- mea.flag2[,c('dsstox_substance_id',
                        'chnm',
                        'casn',
                        'total.endpoints.screened',
                        'up.ct',
                        'dn.ct',
                        'p5.log.acc',
                        'firing',
                        'bursting',
                        'connectivity',
                        'cytotoxicity')]

# should filter based on hitc sum
# should require at least 3-4 to be positive in a single direction
# edit the flag on this basis, then flag becomes the min potency in acute MEA if it is the min potency for the chemical
# consider 5th percentile value of potency after filtering

mea.flag3[, mea.call := 0]
mea.flag3[ up.ct > 3 | dn.ct>3, mea.call := 1]

mea.flag3$min_nam <- asnm.tier1$col_min[match(mea.flag3$dsstox_substance_id,
                                     asnm.tier1$dsstox_substance_id)]
mea.flag3[min_nam=='CCTE', min_nam := 'MEA']
mea.flag3[!min_nam=='MEA', mea.call := 0]

6 Combine potency related flags

6.1 Put MEA and BSK together

head(mea.flag3)

head(bsk.flag2)

mea.bsk.comb <- merge.data.table(mea.flag3[,-c('chnm')],
                                 bsk.flag2,
                                 by=c('dsstox_substance_id'),
                                 fill=TRUE,
                                 all.x=TRUE,
                                 all.y=TRUE)

# deal with issue of different sets of chemicals
# get all numbers into the flag table for potency-related flags

mea.bsk.comb$min_nam <- asnm.tier1$col_min[match(mea.bsk.comb$dsstox_substance_id,
                                     asnm.tier1$dsstox_substance_id)]
mea.bsk.comb[min_nam=='CCTE', min_nam := 'MEA']
mea.bsk.comb$aqc_indicator <- ad.tbl$aqc_indicator[match(mea.bsk.comb$dsstox_substance_id, ad.tbl$DTXSID)]
mea.bsk.comb$apcra.pro.only <- ad.tbl$apcra.pro.only[match(mea.bsk.comb$dsstox_substance_id, ad.tbl$DTXSID)]
mea.bsk.comb$ber.med.aed50 <- asnm.tier1.all.invivo$ber.med.aed50[match(mea.bsk.comb$dsstox_substance_id,
                                                                     asnm.tier1.all.invivo$dsstox_substance_id)]

mea.bsk.comb <- mea.bsk.comb[order(ber.med.aed50)]

6.2 Add HIPPTox organ flags

Here we show the in vitro potency for HIPPTox so it is comparable to MEA and BSK potencies in vitro.

all.organs <- merge.data.table(mea.bsk.comb,
                               asnm.tier1[,c(
                                 'dsstox_substance_id',
                                 'Astar_BEAS2B_final',
                                 'Astar_HepG2_final',
                                 'Astar_HK2_final')],
                               all.x=TRUE,
                               all.y=TRUE,
                               by = c('dsstox_substance_id'))

all.organs$chnm_mea <- mea.ref$chnm[match(all.organs$dsstox_substance_id,
                                              mea.ref$dsstox_substance_id)]
all.organs[is.na(chnm), chnm := chnm_mea]
all.organs[,c('chnm_mea') := NULL]

colnames(all.organs)

6.3 Make target organ hazard flags

In this chunk, we make the table that will go to supplement. Please note that these are all in vitro micromolar potencies, as selective activity and other decisions were made on the in vitro potency scale.

Also note that only 133/196 substances here had MEA data.

#colnames(apcra.list)
nrow(apcra.list[mea.acute.any==1])

## [1] 133

all.organs.flag <- all.organs[,c('dsstox_substance_id',
                                 'chnm',
                                 'apcra.pro.only',
                                 'min_nam',
                                 'ber.med.aed50',
                                 'aqc_indicator',
                                 'p5.log.acc',
                                 'mea.call',
                                 'firing',
                                 'bursting',
                                 'connectivity',
                                 'cytotoxicity',
                                 'acute.tox',
                                 'immun.sup.lec',
                                 'selective',
                                 'Astar_BEAS2B_final',
                                 'Astar_HepG2_final',
                                 'Astar_HK2_final'
                                 )]
setnames(all.organs.flag, c('p5.log.acc',
                                 'mea.call',
                                 'firing',
                                 'bursting',
                                 'connectivity',
                                 'cytotoxicity',
                                 'acute.tox',
                                 'immun.sup.lec',
                                 'selective',
                                 'Astar_BEAS2B_final',
                                 'Astar_HepG2_final',
                                 'Astar_HK2_final'),
         c('MEA_p5_potency',
                                 'MEA_call',
                                 'MEA_firing',
                                 'MEA_bursting',
                                 'MEA_connectivity',
                                 'MEA_cytotoxicity',
                                 'BioMAP_acute',
                                 'BioMAP_immunosupp',
                                 'BioMAP_immuno_sel',
                                 'HIPPTox_Lung',
                                 'HIPPTox_Liver',
                                 'HIPPTox_Kidney')
)

all.organs.flag <- all.organs.flag[order(ber.med.aed50)]
apcra.list[is.na(mea.acute.any), mea.acute.any:=0]
all.organs.flag$mea_avail <- apcra.list$mea.acute.any[match(all.organs.flag$dsstox_substance_id,
                                                            apcra.list$DTXSID)]        

Begin to shape the matrix for heatmap presentation.

organs_mat <- all.organs.flag[apcra.pro.only==1]
organs_mat <- organs_mat[aqc_indicator==1]
organs_mat <- organs_mat[ber.med.aed50 < 4 ]


organs_ref_mat <- unique(all.organs.flag[dsstox_substance_id %in% c(mea.ref.dtxsids, immunosupp.drug.dtxsids)])

chnm_order <- c('Abamectin','beta-Cyfluthrin','Lindane','Tributyltin chloride',
                        'Azathioprine', 'Cyclosporin A', 'Dexamethasone sodium phosphate','Methotrexate')
organs_ref_mat$chnm <- factor(organs_ref_mat$chnm, levels=chnm_order)
organs_ref_mat <- organs_ref_mat[order(as.factor(chnm))]
organs_ref_mat

colnames(organs_mat)

##  [1] "dsstox_substance_id" "chnm"                "apcra.pro.only"     
##  [4] "min_nam"             "ber.med.aed50"       "aqc_indicator"      
##  [7] "MEA_p5_potency"      "MEA_call"            "MEA_firing"         
## [10] "MEA_bursting"        "MEA_connectivity"    "MEA_cytotoxicity"   
## [13] "BioMAP_acute"        "BioMAP_immunosupp"   "BioMAP_immuno_sel"  
## [16] "HIPPTox_Lung"        "HIPPTox_Liver"       "HIPPTox_Kidney"     
## [19] "mea_avail"           "mea_avail_face"      "mea_avail_col"

# comprise the main matrices

# APCRA test chemical matrix
organs_mat2 <- as.matrix(organs_mat[,c("MEA_p5_potency",
                                     #  "MEA_call",
                                       "MEA_firing","MEA_bursting","MEA_connectivity","MEA_cytotoxicity",   
                                       "BioMAP_acute","BioMAP_immunosupp","BioMAP_immuno_sel",
                                       "HIPPTox_Lung","HIPPTox_Liver","HIPPTox_Kidney")])
rownames(organs_mat2) <- organs_mat[, chnm] # define rownames as chemical name

# reference matrix
organs_ref2 <- as.matrix(organs_ref_mat[,c("MEA_p5_potency",
                                           #"MEA_call",
                                           "MEA_firing","MEA_bursting","MEA_connectivity","MEA_cytotoxicity",   
                                       "BioMAP_acute","BioMAP_immunosupp","BioMAP_immuno_sel",
                                       "HIPPTox_Lung","HIPPTox_Liver","HIPPTox_Kidney")])
rownames(organs_ref2) <- organs_ref_mat[, chnm] # define rownames as chemical name

# annotations for main hm
anno_df2 <- data.frame(rownames(organs_mat2))
anno_df2$min_nam <- all.organs.flag$min_nam[match(anno_df2[,1],
                                            all.organs.flag$chnm)]
anno_df2$chnm <- all.organs.flag$chnm[match(anno_df2[,1],
                                      all.organs.flag$chnm)]
anno_df2$ber <- all.organs.flag$ber.med.aed50[match(anno_df2[,1],
                                             all.organs.flag$chnm)]
anno_df2$mea_avail <- all.organs.flag$mea_avail[match(anno_df2[,1],
                                                      all.organs.flag$chnm)]

anno_df2$mea_call <- all.organs.flag$MEA_call[match(anno_df2[,1],
                                                      all.organs.flag$chnm)]

all.organs.flag[mea_avail==1, mea_avail_face := 'plain']
all.organs.flag[mea_avail==0, mea_avail_face := 'italic']
all.organs.flag[MEA_call==1, mea_avail_face := 'bold']

all.organs.flag[mea_avail==1, mea_avail_col := 'black']
all.organs.flag[mea_avail==0, mea_avail_col := 'blue']
all.organs.flag[MEA_call==1, mea_avail_col := 'red']


anno_df2$mea_avail_face <- all.organs.flag$mea_avail_face[match(anno_df2[,1],
                                                      all.organs.flag$chnm)]

anno_df2$mea_avail_col <- all.organs.flag$mea_avail_col[match(anno_df2[,1],
                                                      all.organs.flag$chnm)]

left_ha2 <- rowAnnotation(min_nam = anno_text(anno_df2$min_nam,
                                              gp=gpar(fontface=anno_df2$mea_avail_face,
                                                      col=anno_df2$mea_avail_col,
                                                      fontsize=10,fontfamily='sans'),
                                            just='right',
                                            location=1,
                                            show_name = FALSE))
right_ha2 <- rowAnnotation(BER = anno_barplot(anno_df2$ber, width = unit(4,'cm')))

# annotations for ref hm
anno_ref2 <- data.frame(rownames(organs_ref2))
anno_ref2$chnm <- all.organs.flag$chnm[match(anno_ref2[,1],
                                      all.organs.flag$chnm)]

col_fun = colorRamp2(c(-2, 0, 2, 4, 6), c('#F0F921FF','#E16462FF','#6A00A8FF','#0D0887FF', 'gray'))
hm_organs <- Heatmap(matrix = organs_mat2, 
                       cluster_columns = FALSE,
                       cluster_rows=FALSE,
                       name="log10-mg/kg/day",
                       #col=col_fun,
                       na_col='gray',
                       col = col_fun,
                       show_row_names = TRUE, 
                       row_dend_width = unit(3, "cm"),
                       column_names_max_height = unit(8, "cm"),
                       column_names_gp = gpar(fontsize = 12),
                     row_names_gp = gpar(fontsize=10),
                       heatmap_legend_param = list(title = expression(paste("In vitro potency (",mu,"M)")),
                                                   legend_direction = 'horizontal'),
                       width=unit(6,'cm'),
                       height = unit(14, 'cm'),
                       rect_gp = gpar(col = "gray", lwd = 2),
                       column_names_side='top',
                       left_annotation = left_ha2,
                       right_annotation = right_ha2)

hm_organs_ref <- Heatmap(matrix = organs_ref2, 
                       cluster_columns = FALSE,
                       cluster_rows=FALSE,
                       name="log10-mg/kg/day",
                       col=col_fun,
                       na_col='gray',
                       #col = viridis(option='A', 100),
                       show_row_names = TRUE, 
                       #row_dend_width = unit(3, "cm"),
                       #column_names_max_height = unit(8, "cm"),
                       column_names_gp = gpar(fontsize = 12),
                       row_names_gp = gpar(fontsize=10,fontfamily='sans'),
                       heatmap_legend_param = list(title = expression(paste("In vitro potency (",mu,"M)")),
                                                   legend_direction = 'horizontal'),
                       #width=unit(6,'cm'),
                       height = unit(4, 'cm'),
                       rect_gp = gpar(col = "gray", lwd = 2),
                       column_names_side='top')

hm.org <- draw(hm_organs %v%  hm_organs_ref, 
              heatmap_legend_side='bottom',
              ht_gap = unit(1, "cm"))

7 Make the preliminary output spreadsheet

7.1 Combine flags and AED/in vivo table

This is specific to APCRA chemicals and their data. We are going to add this to an existing table with all columns. Then, we will make a streamlined version of the table for supplement/sharing.

colnames(asnm.tier1.all.invivo)

##  [1] "dsstox_substance_id"                    
##  [2] "aed50.3compss.5pacc"                    
##  [3] "aed50.3compss.atg"                      
##  [4] "aed50.3compss.bsk"                      
##  [5] "aed50.3compss.ccte"                     
##  [6] "aed50.3compss.nvs"                      
##  [7] "aed50.3compss.stm"                      
##  [8] "aed50.3compss.httr.mcf7.sigbmd"         
##  [9] "aed50.3compss.httr.u2os.sigbmd"         
## [10] "aed50.3compss.httr.heparg.sigbmd"       
## [11] "aed50.3compss.htpp.u2os"                
## [12] "aed50.3compss.astar.beas2b"             
## [13] "aed50.3compss.astar.hepg2"              
## [14] "aed50.3compss.astar.hek293"             
## [15] "aed50.pbtk.5pacc"                       
## [16] "aed50.pbtk.atg"                         
## [17] "aed50.pbtk.bsk"                         
## [18] "aed50.pbtk.ccte"                        
## [19] "aed50.pbtk.nvs"                         
## [20] "aed50.pbtk.stm"                         
## [21] "aed50.pbtk.httr.mcf7.sigbmd"            
## [22] "aed50.pbtk.httr.u2os.sigbmd"            
## [23] "aed50.pbtk.httr.heparg.sigbmd"          
## [24] "aed50.pbtk.htpp.u2os"                   
## [25] "aed50.pbtk.astar.beas2b"                
## [26] "aed50.pbtk.astar.hepg2"                 
## [27] "aed50.pbtk.astar.hek293"                
## [28] "the.day"                                
## [29] "aed50.atg"                              
## [30] "aed50.bsk"                              
## [31] "aed50.ccte"                             
## [32] "aed50.nvs"                              
## [33] "aed50.stm"                              
## [34] "aed50.httr.mcf7"                        
## [35] "aed50.httr.u2os"                        
## [36] "aed50.httr.heparg"                      
## [37] "aed50.htpp.u2os"                        
## [38] "aed50.astar.beas2b"                     
## [39] "aed50.astar.hepg2"                      
## [40] "aed50.astar.hek293"                     
## [41] "p5.toxval.numeric"                      
## [42] "p10.toxval.numeric"                     
## [43] "p15.toxval.numeric"                     
## [44] "p20.toxval.numeric"                     
## [45] "p25.toxval.numeric"                     
## [46] "p30.toxval.numeric"                     
## [47] "p5.toxval.numeric.sub"                  
## [48] "p10.toxval.numeric.sub"                 
## [49] "p15.toxval.numeric.sub"                 
## [50] "p20.toxval.numeric.sub"                 
## [51] "p25.toxval.numeric.sub"                 
## [52] "p30.toxval.numeric.sub"                 
## [53] "CASRN"                                  
## [54] "preferred_name"                         
## [55] "apcra.pro.only"                         
## [56] "T0"                                     
## [57] "T4"                                     
## [58] "Call"                                   
## [59] "aqc_iv_pass"                            
## [60] "aqc_indicator"                          
## [61] "AVERAGE_MASS"                           
## [62] "log10VP"                                
## [63] "OCTANOL_WATER_PARTITION_LOGP_OPERA_PRED"
## [64] "logP.indicator"                         
## [65] "mw.indicator"                           
## [66] "logVP.indicator"                        
## [67] "half.life"                              
## [68] "mlr.toxval25.aed50"                     
## [69] "med.aed50"                              
## [70] "min.aed50"                              
## [71] "data_poor_binary"                       
## [72] "ECHA_min_systemic_POD_mkd"              
## [73] "log.ECHA_min_systemic_POD_mkd"          
## [74] "log.repdose.90d.5th"                    
## [75] "log.repdose.any.5th"                    
## [76] "log.repdose.90d.25th"                   
## [77] "log.repdose.any.25th"                   
## [78] "ug.kg.bw.day"                           
## [79] "ttc"                                    
## [80] "chnm"                                   
## [81] "pod.ratio"                              
## [82] "ttc.ratio"                              
## [83] "sub.ratio"                              
## [84] "pod.ratio.size"                         
## [85] "seem3.u95.log10"                        
## [86] "seem3.log10"                            
## [87] "ber.targeted"                           
## [88] "ber.httr"                               
## [89] "ber.htpp"                               
## [90] "ber.astar"                              
## [91] "ber.med.aed50"

apcra.tbl.draft <- merge.data.table(asnm.tier1.all.invivo,
                                    asnm.tier1.all.ratios[,c('dsstox_substance_id',
                                                             'pod.ratio5',
                                                             'pod.ratio25'
                                                             )],
                                    by='dsstox_substance_id',
                                    all.x=TRUE)

Add in critical elements of the DART flag table.

apcra.tbl.draft <- merge.data.table(apcra.tbl.draft,
                                    dart.flag[,c('dsstox_substance_id',
                                                 'chnm',
                                                 'min_nam',
                                                 "dev.test",
                                                 "dev.flag", 
                                                 "dev.flag.specific", 
                                                 "final.er.flag", 
                                                 "final.ar.flag")],
                                    by='dsstox_substance_id',
                                    all.x=TRUE)

Add in critical elements of the target organ system flag table.

apcra.tbl.draft <- merge.data.table(apcra.tbl.draft,
                                    all.organs.flag[,c(
                                      "dsstox_substance_id",
                                      "mea_avail",
                                      "MEA_p5_potency", 
                                    "MEA_call", 
                                    "MEA_firing", 
                                    "MEA_bursting", 
                                    "MEA_connectivity", 
                                    "MEA_cytotoxicity", 
                                    "BioMAP_acute", 
                                    "BioMAP_immunosupp", 
                                    "BioMAP_immuno_sel", 
                                    "HIPPTox_Lung", 
                                    "HIPPTox_Liver", 
                                    "HIPPTox_Kidney")],
                                    by='dsstox_substance_id',
                                    all.x=TRUE)

Calculate the standard deviation of the AED50s by assay source, excluding HIPPTox.

apcra.tbl.draft <- apcra.tbl.draft %>% rowwise() %>% mutate(aed50_sd=sd(c(aed50.atg, aed50.bsk, aed50.ccte, aed50.nvs, aed50.stm, aed50.httr.mcf7,aed50.httr.u2os, aed50.httr.heparg, aed50.htpp.u2os),na.rm=TRUE)) %>% data.table()

#duplicated(apcra.tbl.draft)
apcra.tbl.draft <- unique(apcra.tbl.draft)
#apcra.tbl.draft[c(26:27),]
#apcra.tbl.draft <- apcra.tbl.draft[-c(27),] # remove replicate endosulfan row - give preference to MEA call = 1

Make a reduced table that is easier to navigate.

dput(colnames(apcra.tbl.draft))

apcra.tbl.select <- apcra.tbl.reduced[, c("dsstox_substance_id",
                                        "CASRN", 
                                        "preferred_name", 
                                        "chnm", 
                                        "apcra.pro.only", 
                                        "T0", "T4", "Call", "aqc_iv_pass", "aqc_indicator", 
                                        "AVERAGE_MASS", "log10VP", "OCTANOL_WATER_PARTITION_LOGP_OPERA_PRED", 
                                        "the.day",
                                        "aed50_sd",
                                        "med.aed50",
                                        
                                        "log.p5.toxval.pod", 
                                        "log.p5.toxval.sub",
                                        "log.ECHA_min_systemic_POD_mkd", "log.repdose.90d.pod", "log.repdose.any.pod", 
                                        "ug.kg.bw.day", "ttc", 
                                        "pod.ratio", "ttc.ratio", "sub.ratio", "pod.ratio.size", 
                                        "seem3.u95.log10", 
                                        "seem3.log10", 
                                        "ber.targeted", 
                                        "ber.httr", 
                                        "ber.htpp", 
                                        "ber.astar", 
                                        "ber.med.aed50", 
                                         
                                        "MEA_call", 
                                        
                                        "BioMAP_immunosupp", 
                                        "BioMAP_immuno_sel", 
                                        "HIPPTox_Lung", 
                                        "HIPPTox_Liver", 
                                        "HIPPTox_Kidney") ]

7.2 Save each table

This is to include all reference chemicals associated with each flag

save(ad.tbl,
     #apcra.tbl.reduced,
     #apcra.tbl.select,
     apcra.tbl.draft,
     asnm.tier1.all.ratios,
     dart.flag, 
     er.ar.apcra, 
     dev.flag,
     all.organs.flag,
     bsk.flag2,
     mea.flag3, 
     file='./output/APCRA_haz_flg_BER.RData')

list_data<- list("full_apcra_tbl" = as.data.frame(apcra.tbl.draft),
                 "pod_ratios_apcra_tbl" = as.data.frame(asnm.tier1.all.ratios),
                 "ad.tbl" = as.data.frame(ad.tbl),
          "dart.flag" = as.data.frame(dart.flag),
          "er.ar.apcra" = as.data.frame(er.ar.apcra),
          "dev.flag"= as.data.frame(dev.flag),
          "all.organs.flag" = as.data.frame(all.organs.flag),
          "BioMAP.immune.flag" = as.data.frame(bsk.flag2),
          "MEA.neuro.flag" = as.data.frame(mea.flag3))

write.xlsx(list_data, './output/SuppFile_APCRA_haz_flg_BER_May2024.xlsx')

8 Reproducibility

print(sessionInfo())

## R version 4.2.2 (2022-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 22621)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_United States.utf8 
## [2] LC_CTYPE=English_United States.utf8   
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.utf8    
## 
## attached base packages:
## [1] parallel  grid      stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] ComplexHeatmap_2.14.0 circlize_0.4.15       viridis_0.6.2        
##  [4] viridisLite_0.4.1     tcpl_3.1.0            tidyr_1.3.0          
##  [7] tictoc_1.1            stringr_1.5.1         RMySQL_0.10.25       
## [10] DBI_1.2.2             randomForest_4.7-1.1  plotly_4.10.1        
## [13] openxlsx_4.2.5.2      jtools_2.2.1          kableExtra_1.3.4     
## [16] httk_2.3.0            gridExtra_2.3         gplots_3.1.3         
## [19] ggstance_0.3.6        ggpmisc_0.5.2         ggpp_0.5.2           
## [22] DT_0.28               dplyr_1.1.1           DescTools_0.99.48    
## [25] data.table_1.14.8     cowplot_1.1.1         caret_6.0-94         
## [28] lattice_0.21-8        ggplot2_3.4.2        
## 
## loaded via a namespace (and not attached):
##   [1] readxl_1.4.2         systemfonts_1.0.4    plyr_1.8.8          
##   [4] lazyeval_0.2.2       splines_4.2.2        listenv_0.9.0       
##   [7] digest_0.6.31        foreach_1.5.2        htmltools_0.5.8.1   
##  [10] fansi_1.0.4          memoise_2.0.1        magrittr_2.0.3      
##  [13] cluster_2.1.4        doParallel_1.0.17    confintr_0.2.0      
##  [16] recipes_1.0.5        globals_0.16.2       gower_1.0.1         
##  [19] matrixStats_0.63.0   svglite_2.1.1        hardhat_1.3.0       
##  [22] timechange_0.2.0     colorspace_2.1-0     blob_1.2.4          
##  [25] rvest_1.0.3          mitools_2.4          rbibutils_2.2.13    
##  [28] xfun_0.43            crayon_1.5.2         jsonlite_1.8.4      
##  [31] Exact_3.2            survival_3.5-5       iterators_1.0.14    
##  [34] glue_1.6.2           gtable_0.3.4         ipred_0.9-14        
##  [37] webshot_0.5.4        MatrixModels_0.5-1   GetoptLong_1.0.5    
##  [40] shape_1.4.6          future.apply_1.10.0  BiocGenerics_0.44.0 
##  [43] SparseM_1.81         scales_1.3.0         mvtnorm_1.1-3       
##  [46] Rcpp_1.0.10          tcplfit2_0.1.6       clue_0.3-64         
##  [49] bit_4.0.5            proxy_0.4-27         deSolve_1.35        
##  [52] sqldf_0.4-11         stats4_4.2.2         lava_1.7.2.1        
##  [55] survey_4.1-1         prodlim_2023.03.31   htmlwidgets_1.6.4   
##  [58] httr_1.4.7           RColorBrewer_1.1-3   farver_2.1.1        
##  [61] pkgconfig_2.0.3      nnet_7.3-18          sass_0.4.9          
##  [64] utf8_1.2.3           RMariaDB_1.2.2       labeling_0.4.3      
##  [67] tidyselect_1.2.1     rlang_1.1.0          reshape2_1.4.4      
##  [70] polynom_1.4-1        munsell_0.5.1        cellranger_1.1.0    
##  [73] tools_4.2.2          cachem_1.0.7         cli_3.6.1           
##  [76] gsubfn_0.7           generics_0.1.3       RSQLite_2.3.1       
##  [79] evaluate_0.23        fastmap_1.1.1        yaml_2.3.7          
##  [82] ModelMetrics_1.2.2.2 knitr_1.46           bit64_4.0.5         
##  [85] zip_2.2.2            pander_0.6.5         caTools_1.18.2      
##  [88] purrr_1.0.1          rootSolve_1.8.2.3    future_1.32.0       
##  [91] nlme_3.1-162         quantreg_5.95        xml2_1.3.3          
##  [94] compiler_4.2.2       rstudioapi_0.14      png_0.1-8           
##  [97] e1071_1.7-13         tibble_3.2.1         bslib_0.7.0         
## [100] stringi_1.7.12       highr_0.10           Matrix_1.5-4        
## [103] vctrs_0.6.1          msm_1.7              pillar_1.9.0        
## [106] lifecycle_1.0.4      GlobalOptions_0.1.2  Rdpack_2.4          
## [109] jquerylib_0.1.4      bitops_1.0-7         lmom_2.9            
## [112] R6_2.5.1             KernSmooth_2.23-20   IRanges_2.32.0      
## [115] parallelly_1.35.0    gld_2.6.6            codetools_0.2-19    
## [118] boot_1.3-28.1        MASS_7.3-58.3        gtools_3.9.4        
## [121] chron_2.3-60         proto_1.0.0          rjson_0.2.21        
## [124] withr_3.0.0          S4Vectors_0.36.2     mgcv_1.8-42         
## [127] expm_0.999-7         hms_1.1.3            rpart_4.1.19        
## [130] timeDate_4022.108    class_7.3-21         rmarkdown_2.26      
## [133] pROC_1.18.0          numDeriv_2016.8-1.1  lubridate_1.9.2